
## save stata data file (Rdata file not available at time) to working directory.  Here we preserved the name 2133069031United States LAPOP AmericasBarometer 2017 V1.0_W.dta

#read file
install.packages('readstata13')
library('readstata13')
abus<-read.dta13("2133069031United States LAPOP AmericasBarometer 2017 V1.0_W.dta", convert.factors=F)




library(survey)
library(jtools)
library(weights)



#demographics

abus$sex<-abus$q1-1

abus$ageyears<-abus$q2


abus$hhincome<-(abus$q10new-1)/16

abus$collegedegree[abus$ed_usa %in% c(5, 6)]<-1
abus$collegedegree[!(abus$ed_usa %in% c(5, 6))]<-0



abus$blackvsall[abus$etid==4]<-1
abus$blackvsall[abus$etid!=4|is.na(abus$etid)]<-0

abus$latinovsall[abus$etid==4003]<-1
abus$latinovsall[abus$etid!=4003|is.na(abus$etid)]<-0

abus$asianvsall[abus$etid==4004]<-1
abus$asianvsall[abus$etid!=4004|is.na(abus$etid)]<-0

abus$nativeamvsall[abus$etid==4005]<-1
abus$nativeamvsall[abus$etid!=4005|is.na(abus$etid)]<-0

abus$mixedvsall[abus$etid==4006]<-1
abus$mixedvsall[abus$etid!=4006|is.na(abus$etid)]<-0

abus$othervsall[abus$etid==7|is.na(abus$etid)]<-1
abus$othervsall[abus$etid!=7]<-0

#openness to authoritarian governance items 

#half-sample only
abus$milcoupcrime[abus$jc10==1]<-1
abus$milcoupcrime[abus$jc10==2]<-0

#half-sample only
abus$milcoupcorruption[abus$jc13==1]<-1
abus$milcoupcorruption[abus$jc13==2]<-0


abus$coup1<-abus$milcoupcorruption
abus$coup1[is.na(abus$coup1)]<-0

abus$coup2<-abus$milcoupcrime
abus$coup2[is.na(abus$coup2)]<-0

abus$coup<-abus$coup1+abus$coup2


abus$closecongress[abus$jc15a==1]<-1
abus$closecongress[abus$jc15a==2]<-0

abus$churchillrev<-1-((abus$ing4-1)/6)

abus$disallowvote<-1-(abus$d1-1)/9

abus$disallowdemonstration<-1-(abus$d2-1)/9

abus$disallowrunforoffice<-1-(abus$d3-1)/9

abus$disallowspeeches<-1-(abus$d4-1)/9

abus$toomuchpressfreedombinary[abus$lib1==3]<-1
abus$toomuchpressfreedombinary[abus$lib1 %in% c(1, 2)]<-0

abus$toomuchexpressionfreedombinary[abus$lib2b==3]<-1
abus$toomuchexpressionfreedombinary[abus$lib2b %in% c(1, 2)]<-0

abus$toomuchpolitfreedombinary[abus$lib2c==3]<-1
abus$toomuchpolitfreedombinary[abus$lib2c %in% c(1, 2)]<-0


#economic left items

abus$ownindustry<-(abus$ros1-1)/6

abus$incomeineq<-(abus$ros4-1)/6


#cultural conservatism items

abus$religattend<-1-((abus$q5a-1)/4)

abus$religimport<-1-((abus$q5b-1)/3)

abus$againstsamesexmarriage<-1-((abus$d6-1)/9)


abus$punishcriminals<-(abus$aoj22new-1)/6


#for control variables

abus$trustcongress<-(abus$b13-1)/6

abus$trustparties<-(abus$b21-1)/6

abus$trustlocalgov<-(abus$b32-1)/6

abus$trustelections<-(abus$b47a-1)/6


abus$satisdemoc<-1-((abus$pn4-1)/3)


abus$polint<-1-((abus$pol1-1)/3)

abus$follownews<-1-((abus$gi0-1)/4)





#descriptives


wtd.mean(abus$coup, weights=abus$wt, na.rm=T)

sqrt(wtd.var(abus$coup, weights=abus$wt, na.rm=T))


wtd.mean(abus$closecongress, weights=abus$wt, na.rm=T)

sqrt(wtd.var(abus$closecongress, weights=abus$wt, na.rm=T))


wtd.mean(abus$churchillrev, weights=abus$wt, na.rm=T)

sqrt(wtd.var(abus$churchillrev, weights=abus$wt, na.rm=T))




abus$intol<-with(abus, rowMeans(data.frame(disallowvote, disallowdemonstration, disallowrunforoffice, disallowspeeches)), na.rm=T)

wtd.mean(abus$intol, weights=abus$wt, na.rm=T)

sqrt(wtd.var(abus$intol, weights=abus$wt, na.rm=T))

abusintolitems<-subset(abus, select=c(disallowvote, disallowdemonstration, disallowrunforoffice, disallowspeeches))

psych::alpha(abusintolitems)


abus$toomuchfreedom<-with(abus, rowMeans(data.frame(toomuchpressfreedombinary, toomuchexpressionfreedombinary, toomuchpolitfreedombinary)),na.rm=T)

wtd.mean(abus$toomuchfreedom, weights=abus$wt, na.rm=T)

sqrt(wtd.var(abus$toomuchfreedom, weights=abus$wt, na.rm=T))

abustoomuchfreedomitems<-subset(abus, select=c(toomuchpressfreedombinary, toomuchexpressionfreedombinary, toomuchpolitfreedombinary))

psych::alpha(abustoomuchfreedomitems)


abus$authgov<-  with(abus, rowMeans(data.frame(coup, closecongress, churchillrev, intol, toomuchfreedom)),na.rm=T)


wtd.mean(abus$authgov, weights=abus$wt, na.rm=T)

sqrt(wtd.var(abus$authgov, weights=abus$wt, na.rm=T))

abusauthgovitems<-subset(abus, select=c(coup, closecongress, churchillrev, intol, toomuchfreedom))

psych::alpha(abusauthgovitems)


abus$cultright<-  with(abus, rowMeans(data.frame(religattend, religimport, againstsamesexmarriage, punishcriminals)),na.rm=T)

abus$cultrightnoattend<-with(abus, rowMeans(data.frame(religimport, againstsamesexmarriage, punishcriminals)),na.rm=T)

wtd.mean(abus$cultright, weights=abus$wt, na.rm=T)

sqrt(wtd.var(abus$cultright, weights=abus$wt, na.rm=T))

abuscultrightitems<-subset(abus, select=c(religattend, religimport, againstsamesexmarriage, punishcriminals))

psych::alpha(abuscultrightitems)


abus$econleft<-  with(abus, rowMeans(data.frame(incomeineq, ownindustry)),na.rm=T)

wtd.mean(abus$econleft, weights=abus$wt, na.rm=T)

sqrt(wtd.var(abus$econleft, weights=abus$wt, na.rm=T))

abuseconleftitems<-subset(abus, select=c(incomeineq, ownindustry))

psych::alpha(abuseconleftitems)


abus$confinstit<-with(abus, rowMeans(data.frame(trustelections, trustparties, trustcongress, trustlocalgov)),na.rm=T)

wtd.mean(abus$confinstit, weights=abus$wt, na.rm=T)

sqrt(wtd.var(abus$confinstit, weights=abus$wt, na.rm=T))

abusconfinstititems<-subset(abus, select=c(trustelections, trustparties, trustcongress, trustlocalgov))

psych::alpha(abusconfinstititems)


abus$poleng<-with(abus, rowMeans(data.frame(polint, follownews)),na.rm=T)

wtd.mean(abus$poleng, weights=abus$wt, na.rm=T)

sqrt(wtd.var(abus$poleng, weights=abus$wt, na.rm=T))

abuspolengitems<-subset(abus, select=c(polint, follownews))

psych::alpha(abuspolengitems)



wtd.mean(abus$satisdemoc, weights=abus$wt, na.rm=T)

sqrt(wtd.var(abus$satisdemoc, weights=abus$wt, na.rm=T))



wtd.mean(abus$collegedegree, weights=abus$wt, na.rm=T)

sqrt(wtd.var(abus$collegedegree, weights=abus$wt, na.rm=T))


wtd.mean(abus$ageyears, weights=abus$wt, na.rm=T)

sqrt(wtd.var(abus$ageyears, weights=abus$wt, na.rm=T))



wtd.mean(abus$sex, weights=abus$wt, na.rm=T)

sqrt(wtd.var(abus$sex, weights=abus$wt, na.rm=T))






#main analyses



library(mice)
library(miceadds)
library(survey)

xtract.svymi <- function(z){
  #require(norm)
  m2 <- length(z)
  ns <- nobs(z[[1]]) #first dataset
  betas <- lapply(z, FUN = function(x){coef(x)} )
  sess <- lapply(z, FUN = function(x){vcov(x)}) #the vcov
  #ses <- lapply(sess, FUN = function(x){sqrt(diag(x))}) #the se, not the complete vcov
  
  dm <- lapply(z, FUN = function(x){model.matrix(x)})
  r.1 <- sapply(1:m2, function(row) cor(dm[[row]] %*% betas[[row]], z[[row]]$y))
  fisher.r2z <- function(r) {0.5 * (log(1+r) - log(1-r))}
  zz <- mean(fisher.r2z(r.1))
  r2 <- psych::fisherz2r(zz)^2 #in psych
  R2 <- mean(r2)
  
  test1 <- summary(pool_mi(betas, sess))
  
  test1$term <- row.names(test1)
  test2 <- test1[,c(8,1,2,4)] #4 is the pvalue
  
  row.names(test2) <- c()
  names(test2) <- c('term','estimate','std.error', 'p.value')
  
  tr <- createTexreg(
    coef.names = test2$term, 
    coef = test2$estimate, 
    se = test2$std.error, 
    pvalues = test2$p.value,
    gof.names = c("R2", "Nobs"), 
    gof = c(R2, ns),
    gof.decimal = c(T,F)
  )
  
}

abus<-subset(abus, select=c(sex, ageyears, collegedegree, hhincome, blackvsall, latinovsall, asianvsall, 
                            nativeamvsall, mixedvsall, othervsall, closecongress, coup, churchillrev, disallowvote, 
                            disallowdemonstration, disallowrunforoffice, 
                            disallowspeeches, toomuchpressfreedombinary, toomuchexpressionfreedombinary, 
                            toomuchpolitfreedombinary, 
                            ownindustry, incomeineq, religattend, religimport, againstsamesexmarriage, punishcriminals, 
                            trustelections, trustparties, trustcongress, trustlocalgov, polint, follownews, satisdemoc, estratopri, wt, upm))



abusi<-mice(abus, m=20)

longabus <- complete(abusi, action='long', include=TRUE)

longabus$intol<-with(longabus, rowMeans(data.frame(disallowvote, disallowdemonstration, disallowrunforoffice, disallowspeeches)))

longabus$toomuchfreedom<-with(longabus, rowMeans(data.frame(toomuchpressfreedombinary, toomuchexpressionfreedombinary, toomuchpolitfreedombinary)))

longabus$authgov<-  with(longabus, rowMeans(data.frame(coup, closecongress, churchillrev, intol, toomuchfreedom)))


longabus$cultright<-  with(longabus, rowMeans(data.frame(religattend, religimport, againstsamesexmarriage, punishcriminals)))


longabus$econleft<-  with(longabus, rowMeans(data.frame(incomeineq, ownindustry)))

longabus$age1s<-with(longabus, (ageyears-min(ageyears, na.rm=T))/(max(ageyears, na.rm=T)-min(ageyears, na.rm=T)) )

longabus$confinstit<-with(longabus, rowMeans(data.frame(trustelections, trustparties, trustcongress, trustlocalgov)))

longabus$poleng<-with(longabus, rowMeans(data.frame(polint, follownews)))


longabus$cultrightplussd<-with(longabus, cultright-(mean(cultright, na.rm=T)+sd(cultright, na.rm=T)))

longabus$cultrightminussd<-with(longabus, cultright-(mean(cultright, na.rm=T)-sd(cultright, na.rm=T)))

longabus$econleftplussd<-with(longabus, econleft-(mean(econleft, na.rm=T)+sd(econleft, na.rm=T)))

longabus$econleftminussd<-with(longabus, econleft-(mean(econleft, na.rm=T)-sd(econleft, na.rm=T)))


longabus$cultrightscl <- with(longabus, (cultright-mean(cultright, na.rm=T))/(2*sd(cultright, na.rm=T)))

longabus$econleftscl <- with(longabus, (econleft-mean(econleft, na.rm=T))/(2*sd(econleft, na.rm=T)))

longabus$polengscl <-  with(longabus, (poleng-mean(poleng, na.rm=T))/(2*sd(poleng, na.rm=T)))


longabus$age1sscl <- with(longabus, (age1s-mean(age1s, na.rm=T))/(2*sd(age1s, na.rm=T)))


longabus$satisdemocscl <- with(longabus, (satisdemoc-mean(satisdemoc, na.rm=T))/(2*sd(satisdemoc, na.rm=T)))


longabus$confinstitscl <- with(longabus, (confinstit-mean(confinstit, na.rm=T))/(2*sd(confinstit, na.rm=T)))


longabus$hhincomescl <- with(longabus, (hhincome-mean(hhincome, na.rm=T))/(2*sd(hhincome, na.rm=T)))


longabus$econextrem<-abs(with(longabus, (econleft-mean(econleft, na.rm=T))))

longabus$econextremscl<-with(longabus, (econextrem-mean(econextrem, na.rm=T))/(2*sd(econextrem, na.rm=T)))


abusi2 <- as.mids(longabus)


abus_imp_list <- imputationList(lapply(1:abusi2$m, function(n) mice::complete(abusi2, action = n)))

abusdesign<-svydesign(~upm,strata=~estratopri, nest=TRUE, fpc=NULL, data=abus_imp_list, weights=~wt)


abusreg1<- with(abusdesign,svyglm(authgov~cultrightscl+econleftscl))


abusreg2<- with(abusdesign,svyglm(authgov~cultrightscl+econleftscl+cultrightscl:econleftscl))


abusreg3<- with(abusdesign,svyglm(authgov~cultrightscl+econleftscl+sex+age1sscl+collegedegree+hhincomescl+blackvsall+latinovsall+asianvsall+nativeamvsall+mixedvsall+othervsall+confinstitscl+polengscl+satisdemocscl))


abusreg4<- with(abusdesign,svyglm(authgov~cultrightscl+econleftscl+sex+age1sscl+collegedegree+hhincomescl+blackvsall+latinovsall+asianvsall+nativeamvsall+mixedvsall+othervsall+confinstitscl+polengscl+satisdemocscl+cultrightscl:econleftscl))




abusregright<- with(abusdesign,svyglm(authgov~cultrightplussd+econleftminussd+cultrightplussd:econleftminussd))


abusregleft<- with(abusdesign,svyglm(authgov~cultrightminussd+econleftplussd+cultrightminussd:econleftplussd))


abusregprot<- with(abusdesign,svyglm(authgov~cultrightplussd+econleftplussd+cultrightplussd:econleftplussd))


abusregfree<- with(abusdesign,svyglm(authgov~cultrightminussd+econleftminussd+cultrightminussd:econleftminussd))

abusout1 <- xtract.svymi(abusreg1)
abusout2 <- xtract.svymi(abusreg2)
abusout3 <- xtract.svymi(abusreg3)
abusout4 <- xtract.svymi(abusreg4)


htmlreg(list(abusout1,abusout2,abusout3,abusout4), digits = 3,file = "lapopUS17.html",custom.coef.names = c("Constant","Cultural Conservatism","Left Economic Attitudes","Cultural Conservatism x Left Economic Attitudes","Female","Age","College Degree","Household Income","Black","Latino","Asian","Native American","Mixed Race","Other Race","Confidence in Institutions","Political Engagement","Satisfaction with Democracy"),caption.above = T,caption = "Table D-15: Predictors of Openness to Authoritarian Governance in LAPOP U.S.-2017 Dataset",custom.note = ("%stars. <br> Continuous predictors are mean-centered and scaled by two standard deviations. Binary predictors are coded 0-1.<br>The outcome variable is coded to range from 0.0 to 1.0 <br> Missing data were replaced using multiple imputation with results pooled across 20 copies of the dataset. <br> Comparison Racial-Ethnic Group is White"),custom.gof.names = c("R<sup>2</sup>","Observations"))

abuscultest <- summary(pool(abusreg1),conf.int = T)[2,2]
abuscultlower <- summary(pool(abusreg1),conf.int = T)[2,7]
abuscultupper <- summary(pool(abusreg1),conf.int = T)[2,8]

abuseconest <- summary(pool(abusreg1),conf.int = T)[3,2]
abuseconlower <- summary(pool(abusreg1),conf.int = T)[3,7]
abuseconupper <- summary(pool(abusreg1),conf.int = T)[3,8]

abusintest <- summary(pool(abusreg2),conf.int = T)[4,2]
abusintlower <- summary(pool(abusreg2),conf.int = T)[4,7]
abusintupper <- summary(pool(abusreg2),conf.int = T)[4,8]


abusrightest <- summary(pool(abusregright),conf.int = T)[1,2]
abusrightlower <- summary(pool(abusregright),conf.int = T)[1,7]
abusrightupper <- summary(pool(abusregright),conf.int = T)[1,8]


abusleftest <- summary(pool(abusregleft),conf.int = T)[1,2]
abusleftlower <- summary(pool(abusregleft),conf.int = T)[1,7]
abusleftupper <- summary(pool(abusregleft),conf.int = T)[1,8]


abusprotest <- summary(pool(abusregprot),conf.int = T)[1,2]
abusprotlower <- summary(pool(abusregprot),conf.int = T)[1,7]
abusprotupper <- summary(pool(abusregprot),conf.int = T)[1,8]


abusfreeest <- summary(pool(abusregfree),conf.int = T)[1,2]
abusfreelower <- summary(pool(abusregfree),conf.int = T)[1,7]
abusfreeupper <- summary(pool(abusregfree),conf.int = T)[1,8]


abuscultb3<-extract.cis(abusreg3)[2,1]
abuspolengb3<-extract.cis(abusreg3)[15,1]
abusageb3<-extract.cis(abusreg3)[5,1]
abuscollegeb3<-extract.cis(abusreg3)[6,1]
abusconfinstitb3<-extract.cis(abusreg3)[14,1]



#additional analyses without imputation, scales with all available items and listwise deletion in analysis. Reported in Part E of SOM.

abus$cultrightscl<-with(abus, (cultright-mean(abus$cultright, na.rm=T))/(2*sd(cultright, na.rm=T)))

abus$econleftscl<-with(abus, (econleft-mean(abus$econleft, na.rm=T))/(2*sd(econleft, na.rm=T)))


abus$cultrightplussd<-with(abus, cultright-(mean(cultright, na.rm=T)+sd(cultright, na.rm=T)))

abus$cultrightminussd<-with(abus, cultright-(mean(cultright, na.rm=T)-sd(cultright, na.rm=T)))

abus$econleftplussd<-with(abus, econleft-(mean(econleft, na.rm=T)+sd(econleft, na.rm=T)))

abus$econleftminussd<-with(abus, econleft-(mean(econleft, na.rm=T)-sd(econleft, na.rm=T)))


abusdesign<-svydesign(~upm, strata=abus$estratopri, nest=TRUE, fpc=NULL, data=abus, weights=abus$wt)


reg1abusl<- svyglm(authgov~cultrightscl+econleftscl, design=abusdesign)


reg2abusl<- svyglm(authgov~cultrightscl+econleftscl+cultrightscl:econleftscl, design=abusdesign)

abuslregright<- svyglm(authgov~cultrightplussd+econleftminussd+cultrightplussd:econleftminussd, design=abusdesign)


abuslregleft<- svyglm(authgov~cultrightminussd+econleftplussd+cultrightminussd:econleftplussd, design=abusdesign)


abuslregprot<- svyglm(authgov~cultrightplussd+econleftplussd+cultrightplussd:econleftplussd, design=abusdesign)


abuslregfree<- svyglm(authgov~cultrightminussd+econleftminussd+cultrightminussd:econleftminussd, design=abusdesign)



abuslcultest<-summary(reg1abusl)$coef[2,1]
abuslcultlower<-confint(reg1abusl)[2,1]
abuslcultupper<-confint(reg1abusl)[2,2]

abusleconest<-summary(reg1abusl)$coef[3,1]
abusleconlower<-confint(reg1abusl)[3,1]
abusleconupper<-confint(reg1abusl)[3,2]

abuslintest<-summary(reg2abusl)$coef[4,1]
abuslintlower<-confint(reg2abusl)[4,1]
abuslintupper<-confint(reg2abusl)[4,2]


abuslrightest<-summary(abuslregright)$coef[1,1]
abuslrightlower<-confint(abuslregright)[1,1]
abuslrightupper<-confint(abuslregright)[1,2]

abuslleftest<-summary(abuslregleft)$coef[1,1]
abuslleftlower<-confint(abuslregleft)[1,1]
abuslleftupper<-confint(abuslregleft)[1,2]

abuslprotest<-summary(abuslregprot)$coef[1,1]
abuslprotlower<-confint(abuslregprot)[1,1]
abuslprotupper<-confint(abuslregprot)[1,2]

abuslfreeest<-summary(abuslregfree)$coef[1,1]
abuslfreelower<-confint(abuslregfree)[1,1]
abuslfreeupper<-confint(abuslregfree)[1,2]

